arXiv:1507.06795vl [cond-mat.soft] 24Jul2015 


Self-assembly of colloidal bands driven by a periodic external field 

Andre S. Nunes/ Nuno A. M. Araujo/’[^ and Margarida M. Telo da Gama^ 

^ Departamento de Fisica, Faculdade de Ciencias, Universidade de Lisboa, 

P-1749-016 Lisboa, Portugal, and Centro de Fisiea Teoriea e Computaeional, 
Universidade de Lisboa, P-1749-016 Lisboa, Portugal 

We study the formation of bands of colloidal particles driven by periodic external fields. Using 
Brownian dynamics, we determine the dependence of the band width on the strength of the particle 
interactions and on the intensity and periodicity of the held. We also investigate the switching 
(held-on) dynamics and the relaxation times as a function of the system parameters. The observed 
scaling relations were analyzed using a simple dynamic density-functional theory of huids. 


INTRODUCTION 

The possibility of obtaining materials with enhanced 
physical properties from the self-assembly of colloidal 
particles has motivated experimental and theoretical 
studies over decades mM- Initially, the focus was on 
identifying novel phases and constructing the equilibrium 
phase diagrams, based on the properties of the individual 
particles (shape, size, and chemistry) [IHZ]- However, the 
impressive advance of optical and lithographic techniques 
opened the possibility of exploring alternative routes as, 
for example, the use of substrates and interfaces IHHI3], 
the control of the suspending medium m or the assem¬ 
bly under flow More recently, there has been a sus¬ 
tained interest on the self-assembly under the presence 
of an electromagnetic (EM) field p^I[2Q] . Uniform EM 
fields couple to the rotational degrees of freedom of mag¬ 
netic aspherical particles allowing to control their orien¬ 
tation [21] or fine-tune the strength and directionality of 
the particle interactions ng HOI Eg. Space-varying pe¬ 
riodic fields are employed to impose constraints on the 
particle position, forming virtual molds that induce spa¬ 
tial periodic patterns [23H25] . 

Under external constraints, colloidal suspensions are 
usually driven out of equilibrium and the relaxation to¬ 
wards equilibrium results from the competition of various 
mechanisms occurring at different length and time scales. 
Thus, besides the identification of the equilibrium struc¬ 
tures and their dependence on the experimental condi¬ 
tions, it is of paramount practical interest to character¬ 
ize the kinetic pathways towards the desired structures 
and the timescales involved. Here, we study the proto¬ 
typical example of field-driven self-assembly of colloidal 
particles. Eor simplicity, we consider the formation of col¬ 
loidal bands driven by a periodic (sinusoidal) electromag¬ 
netic field and study how the equilibrium structures and 
the dynamics towards equilibrium depend on the field 
strength and periodicity, as well as on the particle inter¬ 
actions. We perform extensive Brownian dynamics simu¬ 
lations and complement the study with a coarse-grained 
analysis based on the dynamic density-functional theory 
of fluids (DDET). 

The paper is organized in the following way. In Sec¬ 


tion [] we present the model and simulation details. The 
dependence of the stationary state structure and relax¬ 
ation dynamics on the model parameters is discussed in 
Section [] Einally, we draw some conclusions in Section [] 


MODEL AND SIMULATIONS 

We consider a two dimensional system of colloidal par¬ 
ticles in the overdamped regime. The pairwise particle 
interactions are described by the repulsive Yukawa po¬ 
tential. 
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where r = \ri — rj\ is the distance between particles i 
and j. Vb sets the energy scale and the screening param¬ 
eter a sets the range of the potential. A characteristic 
particle radius, r^, can be defined as Vp = {2a)~^, which 
we set as the unit of length. Eor simplicity, we consider 
only repulsive interactions (Vb >0). To quantify the 
pairwise interactions we define the integrated parameter 
A = f Vij{r)dr expressed in units of ksTUp. In the nu¬ 
merical simulations we set a = 1/2 and change A by 
varying Vb- 

We consider an external field that is constant along 
the ^-direction and periodic along the x-direction. The 
particle/fleld interaction is described by a periodic po¬ 
tential. 


Vext{x, y) = Vext{x) = Ve sin {nx ), 
where Ve is the strength of the potential and 
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is the wave number. P sets the number of minima of the 
potential in a system of linear length L. 

We perform Brownian dynamics simulations where the 
equation of motion of the colloidal particle i is 
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FIG. 1. Snapshot of the system for t = 0, t = 4 and t = 
16. The plots correspond to the density profiles along the 
x-direction. Under the influence of the external potential the 
particles organize into bands. The linear length of the system 
is L — 100, the interaction parameters are Vb = 5 and a — 
1/2. The strength of the field is Ve = 5 and the number of 
minima is P = 4. 

where 7 is the Stokes friction coefficient and the last term 
on the right-hand side ^i{t) is the Langevin force mim¬ 
icking the particle/fluid interaction, which is randomly 
sampled from a Gaussian distribution of zero mean and 
second moment = 2kBTj6kiS{t — t'). The 

indices k and I run over the spatial dimensions, ks is 
the Boltzmann constant and T the thermostat tempera¬ 
ture of the suspending fluid. Thus, the random force is 
uncorrelated in time and space. 

We start our simulations with a uniform (random) 
distribution of colloidal particles and switch on the 
field at t = 0. The simulations were performed in¬ 
side square boxes of three different linear lengths L = 
{100,120,150}, in units of the particle radius. By 
fixing the density p = 0.0696, we simulate N = 
{696,1000,1560} colloidal particles, respectively. 

Hereafter, the strength of the external potential Ve is 
expressed in units of ksT and time is defined in units of 
the Brownian time r = r‘^^{kBT)~^^ which is the time 
over which a colloidal particle diffuses over a region equiv¬ 
alent to its area. To integrate the stochastic differential 
equations of motion we followed the scheme proposed 
by Brahka and Heves[26]. which consists of a second- 
order stochastic Runge-Kutta scheme, with a time-step 
of At = 10 “^r. 


RESULTS 

In the absence of external fields the spatial distribu¬ 
tion of the colloidal particles is homogeneous. When the 
external field is switched on the particles are dragged 
towards the closest minimum and self-organize into in¬ 
dividual bands around it (see snapshots in Fig. §. To 


quantify the spatial arrangement of the colloidal parti¬ 
cles we measure the density, p(r), defined as the number 
of particles per unit area. Given the symmetry of the 
external field (see Eq. we focus on the profile of the 
density along the x-direction. Figure shows a snapshot 
and the density profile for a system with L = 100 at three 
different times. The density profile has the symmetry of 
the external potential. The maxima of the density cor¬ 
respond to the minima of the potential. Similarly, the 
density vanishes at the maxima of the potential. We run 
the simulations for 6 x 10^ timesteps. After this time, 
the changes in the density profile are within the error 
bars and thus we consider that the stationary state has 
been reached. Simulations with different system sizes 
were preformed and we found no significant finite-size 
effects. 

The collective dynamics and the band structure re¬ 
sult from a complex interplay of the particle/fluid, par- 
ticle/field and particle/particle interactions. Here, we 
analyze how the final structure and dynamics depend on 
the model parameters. We start with an analysis of the 
stationary state and proceed with the study of the dy¬ 
namics. 

Stationary State 

Numerical simulations. 

We first characterize the stationary state that corre¬ 
sponds to the thermodynamic equilibrium configurations. 
For P > Ll2 the average particle/field interaction van¬ 
ishes, as the wave length of the potential is shorter than 
the particle diameter. The limit of no external field is 
then recovered. For P < I//2 and strong enough field, 
particles accumulate around the minima forming a band 
structure. Thermal noise and particle/particle interac¬ 
tion cause a broadening of the bands, giving them an 
effective thickness. This thickness is a property of the 
equilibrium configuration. To quantify it, we measure 
the mean square displacement of particles around the lo¬ 
cal minimum 

O' = {xl) - {Xbf (5) 

where the average is over particles in the same band, 
i.e., in between two consecutive maxima of the potential. 
Xl) is their position relative to the local minimum along 
the X-direct ion. The average displacement is zero as the 
distribution of particles is symmetrical with respect to 
the center of the band (minimum of the potential). 

Figure (left panel) shows the dependence of a on the 
model parameters. One clearly sees that, for a given ther¬ 
mostat temperature, the thickness of the bands results 
from the competition between particle/field and parti¬ 
cle/particle interactions. The particle/field interaction 
promotes the formation of thin bands. The stronger the 
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FIG. 2. Mean square displacement of the particle position around the potential minima. On the left-hand side are the results 
obtained from the Brownian dynamics simulations and on the right-hand side those from the DFT calculations. The data is 
obtained by averaging over 25 samples rescaled by P~^. Black, blue and red represent respectively Ve = {3,5,7} and the 
circle, square, triangle and diamond symbols correspond respectively to P = {1,2, 3,4}. The error bars are smaller than the 
size of the symbols. The particle/particle and particle/field interactions have opposite effects on the final configuration. While 
the former promotes wider bands the latter favors thinner ones. 


field (Ve) the thinner the band is. By contrast, the par¬ 
ticle/particle repulsion tends to homogenize the spatial 
distribution of the particles and favors wider bands. Con¬ 
sequently, a increases monotonically with A. Note that 
when we rescale a by P^ a data collapse is obtained for 
different P values into a single curve that depends only 
on Ve and A. In the next section we perform a density- 
functional theory (DFT) analysis to study these depen¬ 
dences. 


Density-functional theory analysis. 

An approximation for the Helmholtz free energy func¬ 
tional can be written as 

r[p\ = J kBTp{r) [log [p{f)h?) - l] df 

+ fy J p{r)p{r)Vpp{r-f*)df'df (6) 
+ J p{^Vext{r)df. 

On the right-hand side, the first term is the free en¬ 
ergy of the ideal gas where A is the thermal de Broglie 
wavelength, the second term is the mean-field approx¬ 
imation to the contribution from the interactions and 
the last term is the external potential contribution. We 
use the local density approximation (LDA) by setting 
p{r^) ^ p{r) which is a good assumption if the density is 
a smooth function of the position. 

Now we employ the definition of chemical potential, 



FIG. 3. Density prohle for L = 100, Ve = 5, A = 43.71 
and P = 4. The dots are the numerical results obtained by 
averaging over 100 samples and the line is the result of DFT. 
The theory is in good agreement with the simulations but it 
underestimates the density at the maxima. 

^ obtain 

kBT\og{p{f)A‘^) + Ap{r) +Vext(,r) = p. (7) 

where A is the interaction parameter defined previously. 
Rearranging the last equation we find an expression for 
the density profile 

p{x) = Z exp[—/3 (Ve sin(ra) + Ap{x))] (8) 

where ^ and Z = A“^e^^ is a normalization 

constant that represents the density of an ideal gas with 
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chemical potential {i. Recall that from the symmetry 
of the external potential the density is only a function 
of the x-coordinate. To solve this equation we impose 
an additional constraint N = J p{r)df arising from the 
conservation of the total number of particles. 



FIG. 4. Maximum density of the profile. Each data point is 
an average over 25 samples and the lines are the DFT results. 
Black, blue and red represent respectively Ve — {3, 5, 7} and 
the circle, square, triangle and diamond symbols correspond 
respectively to P = {1,2, 3,4}. The error bars are smaller 
than the size of the symbols. Note that the density does not 
depend on P. The deviation of the theoretical curves from 
the simulation results increases particle interactions increase. 
The effect of A and Ve on the maximal density is opposite to 
their effect on the mean square displacement. 


Figure shows the density profile obtained from nu¬ 
merical simulation and DFT for the same set of parame¬ 
ters. The DFT results deviate from the simulation only 
at the maximal and minimal densitities. The maximal 
density is underestimated and the minimal is slightly 
overestimated. We evaluated how this deviation depends 
on the model parameters. Figure|^shows the dependence 
on A of the maximal density pmax ^ defined as the density 
in the center of the band, for different values of Ve and P. 
We find that, while the deviation of the DFT calculation 
from the simulation increases slightly with A it does not 
vary significantly with Ve and P. This also leads to the¬ 
oretical values of the mean square displacement cr, which 
are larger than those obtained by simulation. However, 
as shown in Fig. (right panel) one recovers the same 
qualitative dependence on the model parameters and the 
values of a are of the same order of magnitude. 

Since the density profile is symmetric with respect to 
the center of the band, the mean square displacement is 
given by. 


a = 



x^p''{x)dx. 


(9) 


where p*(x) = p{x)N^ ^ is the probability density func- 

L 

tion for particles in the band and = L p(x)dx. 

2P 

The L outside of the integral accounts for the integration 
of the domain along the ^-direction. Assuming that p(x) 
is constant and Nf) = N/P^ 


a = 


P2 

12P2 


( 10 ) 


Thus, (7 ^ P ^ as observed both from the simulations 
and the DFT analysis. Note that the maximal density (in 
Fig.g does not depend on P. From Eq. we also con¬ 
clude that, near the minima, the external potential and 
the interaction terms have different signs. Consequently, 
Ve and A have opposite effects on a, as observed in the 
simulations and DFT calculations, see Fig. 

By taking the logarithm of Eq. @ we obtain 


log(p(x)) -h l3Ap{x) = log{Z) — pVEsin{nx). (11) 


The dependence of the first term of the left-hand side on 
the density can be neglected with respect to the second 
term in two limiting cases: for strong particle/particle 
interactions or high enough densities. In the former, the 
strong particle/particle repulsion hinders the formation 
of bands and homogenizes the density. The external po¬ 
tential acts as a perturbation to the uniform distribution 
of particles, which is sinusoidal in space and linear on 
Ve‘ The contributions from the external potential and 
particle interactions are given by p{x) which means that 
the dependence of the mean square displavement on the 
strength of the potential intensity is the same in both 
limits, as seen in Fig. [^a) and (b). For weak particle in¬ 
teractions the curves are linear when the strength of the 
potential is also small. But as the particle/particle in¬ 
teraction increases, the potential intensity where the de¬ 
pendence deviates from linear also increases. Note that 
the curve for A = 0 corresponds to the ideal gas where 
the dependence of p{x) on the strength of the potential 
is exponential as expected from Eq. For = 0 the 
density is uniform and the value of a is given by Eq. (10) 
and is the same for all curves. 


Relaxation Dynamics 

Numerical simulations. 

We analyze now the relaxation towards the stationary 
state shown in Fig. Namely, we investigate how the 
relaxation time depends on the model parameters. 

We start with a uniform spatial distribution of par¬ 
ticles, p{x) = po- As the external potential is switched 
on the colloidal particles move accordingly and regions of 
high density of colloids are formed (bands). For a system¬ 
atic analysis, we focus on the evolution of the positions of 
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FIG. 5. DFT results for the mean square displacement for 
P = 2. The lines in (a) are for A = {0,20,40,60,80,100} 
and the lines in (b) are for different initial densities po = 
{2.50,6.96,10.0,15.0,20.0,25.0} x 10“^. For strong inter¬ 
particle interactions and weak external potentials the depen¬ 
dence is linear in Ve. The linear dependence is also observed 
at high densities. 


the minima, Pmin^ corresponding to the potential max¬ 
ima, and the positions of maxima, pmax^ corresponding 
to the potential minima. Figure shows the convergence 
of these densities to their value in the stationary state, 
Poo, foi* a particular set of parameters. One clearly sees 
an exponential decay in time. 

We define the relaxation time as the inverse of 

the slopes of the evolution curves. Figures and show 
that the dependence on the model parameters is differ¬ 
ent for the minima and the maxima suggesting that the 
underlying relaxation mechanisms are different. In gen¬ 
eral, strong potentials favor short relaxation times while 
the dependence on the particle/particle interaction is not 
straightforward. Also, data collapse is observed when we 
rescale r with P^. We will now use DDFT to shed light 
on these findings. 


Dynamic density-functional theory analysis. 

If we consider that the system evolves adiabatically, 
the evolution equation of the local density can be writ- 



FIG. 6. Evolution of the maximal and minimal densities. 
The relaxation into the stationary state is exponential. The 
scattered data are the numerical results and the lines the ex¬ 
ponential fitting. The results shown are averages over 100. 
The parameters are L = 100, Ve = 5, A = 43.71 and P = 4. 


ten directly from the equilibrium Helmholtz free energy 
functional [28] . 


^ dt 


p{f,t)W 


5r[p{f, t)] 
5p{f, t) 


( 12 ) 


Using the functional defined in Eq. (|^, we obtain a 
diffusion equation for the density, 


dp 


= V. I^ApVp + Vei^cos{kx)p 


, (13) 


where the first and last terms are related to the parti¬ 
cle/particle and particle/fiuid interactions, respectively. 
They tend to smooth out any spatial variation of the den¬ 
sity. The middle term is related to the particle/field inter¬ 
action and promotes the formation of bands. Note that, 
when one applies V. to the first term one gets ApV^p, a 
non-linear diffusion term where the coefficient increases 
with the density. It is the interplay of the thermostat 
temperature, external potential and particle/particle in¬ 
teractions that controls the kinetics of relaxation and de¬ 
fines the relaxation time towards the stationary state. A 
similar diffusion equation can also be obtained from the 
Fokker-Planck formalism [29] . 

We focus on the maximal and minimal densities where 
Vp = 0. Let us assume that = p 2 is constant in 
time at these points. At the minimum p 2 > 0 and at the 


maximum p 2 < 0. By solving Eq. (13) we get 


Pmin (^) — 


kBTp2 
Ap2 T Vei^‘^ 


+ Ki exp 


Ap2 T Ve^‘'" 
7 


(14) 

where the expression with the minus sign is for the den¬ 
sity at the minimum and the plus sign for the density at 
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FIG. 7. Relaxation time of the minimal density. The scat¬ 
tered data are the simulation results and the solid lines the 
ones from the DDFT calculations. Black, blue and red rep¬ 
resent respectively Ve = {3,5,7} and the circle, square, tri¬ 
angle and diamond symbols correspond respectively to P = 
{1,2, 3,4}. The simulation results are averages over 45 sam¬ 
ples for Ve — and 25 samples in the other cases. The relax¬ 
ation time does not vary monotonically with the interaction 
strength. It exhibits a maximum that depends on the bal¬ 
ance between A and Ve- The relaxation time is proportional 
to 



FIG. 8. Relaxation time of the maximal density. The scat¬ 
tered data are the numerical results and the solid lines the 
ones from the DDFT calculations. Black, blue and red rep¬ 
resent respectively Ve = {3, 5, 7} and the circle, square, tri¬ 
angle and diamond symbols correspond respectively to P = 
{1,2, 3,4}. The simulation results are averages over 45 sam¬ 
ples for Ve — and 25 samples in the other cases. For strong 
interparticle interactions the relaxation time depends weakly 
on the external potential. 


the maximum. By fixing the initial density it is possible 
to calculate the constant 


K = A 4. ^bTp2 
^ ~ L2 Ap2 =F ■ 

If we define the characteristic relaxation time as 

p{t) - Poo 


(15) 


(16) 


where poo is the density at the given point at equilibrium, 
we arrive at the following expression 


relax 

'min 

max 


7 


Ap2 tVek^' 


(17) 


In the region of parameters where = const the re- 
laxation time is inversely proportional to A and Ve- The 
friction coefficient also plays a role on the dynamics and 
the relaxation time is a monotonic increasing function 
of the viscosity. Replacing k using Eq. ^ we conclude 
that P“^. The relaxation time decreases with 

the number of bands as, on average, particles have to 
travel a shorter distance to reach the nearest potential 


minimum. Eq. (13) is solved numerically using first or¬ 


der finite elements in the 2D square domain with peri¬ 
odic boundary conditions [30]. Eigures[^ and show the 
relaxation time of the minimal and maximal densities, 
respectively. They collapse for P~^ but the inverse pro¬ 
portionality is only observed for high A which is where 
our approximation holds. 


The theoretical results underestimate the relaxation 
time by comparison with the data from the Brownian dy¬ 
namics simulations, consistent with the fact that DDET 
overestimates the diffusion coefficient m- The devia¬ 
tion is greater at low Ve and high A, when the parti¬ 
cle/particle interactions dominate the dynamics. Recall 
that the term for particle/particle interaction in Eq. §) 
is the only approximate term. 

Note that the relaxation time is always positive. Eor 
the minimal density the curvature p 2 is typically small 
and the external potential term dominates. Eor the max¬ 
imal density the curvature is negative and the parti¬ 
cle/particle interaction term dominates over the external 
potential. In the minimum the particles try to escape 
and so the external potential is the relevant interaction 
leading to a decrease in the density. However, the process 
is still influenced by the presence of neighboring particles 
and the relaxation time increases with A at low Eig. 

Eor strong interparticle repulsions there is a broad¬ 
ening of the bands; the particles will then travel shorter 
distances to reach their equilibrium positions, decreasing 
r'^eiax minimal density. By contrast, for the max¬ 

imal density the repulsions strongly affect the dynamics. 

The relaxation times obtained from the simulations 
and the theory coincide for strong potentials as the ef¬ 
fect of the interparticle interactions is irrelevant. When 
Ve decreases the particle/particle repulsion becomes rel¬ 
evant and the EDA approximation leads to increasing 
deviations from the simulation results. Equation ( [T7| ) 
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suggests that for high enough values of A the relaxation 
times do not depend on Ve and P and converge to the 
same value in line with the results shown in Figs, [^and 
For lower A at the maxima, where the particles are 
compressed, equilibrium takes longer to be reached. 

CONCLUSIONS 

We studied the dependence of the structure and re¬ 
laxation dynamics of the equilibrium state of a colloidal 
suspension on the external potential and the interpaticle 
interactions. The analysis based on density-functional 
theory was found to be in good agreement with the re¬ 
sults from Brownian dynamics simulations. 

In the stationary state the system exhibits a band-like 
structure. The thickness of the bands measured in terms 
of the mean square displacement, cr, decreases with the 
strength of the potential, which leads in turn to higher 
densities around the external potential minima. We also 
showed that a scales as where P is the number 

of bands. In the limit of weak external potentials, a 
increases linearly with the potential strength, Ve- 

The relaxation towards the stationary extreme densi¬ 
ties is exponential but the relaxation times are different. 
The system takes longer to reach equilibrium at the den¬ 
sity maxima except when the external potential becomes 
significant weak in which case the relaxation time does 
not depend on the position. At the density maxima 
decreases when both Ve and A increase but at the density 
minima it displays a more complex behavior and there is 
a maximum relaxation time that depends on the balance 
between Ve and A. 

Electromagnetic fields are used in the laboratory to 
promote the formation of a variety of colloidal structures. 
Our results quantify the relaxation times and the equilib¬ 
rium structures of simple colloidal systems, which exhibit 
a rich behaviour that may be controlled in a straightfor¬ 
ward manner. 
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